function y=Hessian(beta, Data)

N=max(Data(:,1));
task=max(Data(:,2));
option=max(Data(:,3));
L=size(beta,2);
y=zeros(L,L);


for id=1:N
   for taskid=1:task
      if sum(Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option,4))~=12
        options=Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option,3);
        rank=Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option,4);
        w=Data((id-1)*task*option+(taskid-1)*option+1:(id-1)*task*option+taskid*option-1,5:5+L-option);
        w=[eye(option-1) w]; 
        
        %define the numerators of the likelihood
        numer(1:option-1,1)=exp(w*beta');

        %define the denominator 
        denom=sum(numer)+1;
        
        %set up the likelihood vector
        probi=numer/denom;
        
        %set up the dprob/dbeta vector (option-1 by L matrix)
        
        %set up a vector to be used in the SOC
        for l=1:L
            expprobi(l)=w(:,l)'*numer/denom;
        end
        for op=1:option-1
           dpdb(op,:)=probi(op)*(w(op,:)-expprobi); 
        end
        
        %set up the SOCs
        y=y+w'*dpdb;
      end
   end
end

y=inv(y);